Variational QMC study of a Hydrogen atom in jellium 
with comparison to LSDA and LSDA-SIC solutions 
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A Hydrogen atom immersed in a finite jellium sphere is solved using variational quantum Monte 
Carlo (VQMC). The same system is also solved using density functional theory (DFT), in both 
the local spin density (LSDA) and self-interaction correction (SIC) approximations. The immersion 
energies calculated using these methods, as functions of the background density of the jellium, are 
found to lie within leV of each other with minima in approximately the same positions. The DFT 
results show overbinding relative to the VQMC result. The immersion energies also suggest an 
improved performance of the SIC over the LSDA relative to the VQMC results. The atom-induced 
density is also calculated and shows a difference between the methods, with a more extended Friedel 
oscillation in the case of the VQMC result. 

PACS numbers; 02.70.Ss, 71.10.Ca, 71.15.Mb, 31.15.Ew 
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I. INTRODUCTION 

The model system of an atom embedded in infinite jel- 
lium is of interest because it can be used to describe the 
bonding of an atom within a solid. The idea is that the 
atom is embedded in a homogeneous electron gas which 
is regarded as the sum of the density tails from all other 
atoms in the solid spatially averaged over the unit cell 
of the atom in question. The immersion energy of the 
atom (the total energy of the atom in jellium minus the 
energies of the separate atom and jellium) can then be 
regarded as a first order approximation to the cohesive 
energy per atom of the solid. In fact, within the frame- 
work of the effective medium theory (EMT)fii^ one in- 
cludes an additional term in this cohesive energy which 
describes the electrostatic attraction between the density 
tails within the unit cell and the Hartree potential set up 
by the atom. This allows one to calculate cohesive prop- 
erties of solids across the periodic table and indeed one 
finds that the Wigner Seitz radii, bulk moduli and cohe- 
sive energies follow the same trends as the experimental 
values. 

In this paper the atom in jellium system is solved 
using the variational quantum Monte Carlo (VQMC) 
method^i^ and the local spin density (LSDA)^i^iii^i^ and 
self- interaction correction (SIC)^" approximations of den- 
sity functional theory (DFT). The particular system used 
is a Hydrogen atom immersed in jellium. In order to solve 
this system using VQMC, the system must be of finite 
size and therefore we approximate the atom in infinite 
jellium with an atom in a finite jellium sphere centered 
on the atom. The quantities to be calculated include 
the total energies of both the Hydrogen in jellium and 
the jellium sphere by itself, the immersion energy and 
the atom-induced density (the density of the Hydrogen 
atom in jellium minus the density of the jellium). These 
quantities are to be calculated as a function of the posi- 
tive background density of the jellium, no, for a jellium 
sphere with a fixed number of electrons. 



VQMC calculations have already been attempted on 
the system of a Hydrogen atom in jellium by Sugiyama 
et al.^^ These calculations however resulted in immer- 
sion energies which differed significantly from the LDA 
results. Atom-induced densities calculated by Sugiyama 
et al. were also found to differ significantly between these 
methods. Our results are based upon a careful study of 
size effects and show a much closer agreement between 
LSDA and QMC. 

The system of a jellium sphere with no embedded atom 
has been solved using QMC (in both the variational and 
diffusion variants) by Sottile et alM- Jellium spheres with 
up to 106 electrons were reported and good agreement 
was found between the LSDA and the diffusion QMC 
(DQMC) results. 

Work on using a Hydrogen atom in a finite jellium 
sphere has also been reportedi^ within the LSDA. It was 
shown that by controlling the size of the jellium sphere, 
a good approximation to Hydrogen in infinite jellium can 
be obtained. 

SIC calculations have already been reported for the 
atom in infinite jellium system.-'^- The EMT was applied, 
and for atoms up to and including the 2p elements the 
SIC cohesive energy was found to be higher than that 
of the LSDA. The interpretation placed on this was that 
the SIC was correcting for the overbinding present in 
the LSDA. Here we compare SIC with both LSDA and 
VQMC. 

In the following section we first investigate the LSDA 
solution of a Hydrogen atom in a finite jellium sphere. 
We examine the effect of changing the jellium sphere size 
and establish choices of sphere which yield good approxi- 
mations to the infinite jellium solution. We then consider 
in Sections |TTT] and |IV] the VQMC and SIC methods for 
the same finite sphere system. Section |V] presents a com- 
parison of our VQMC results with those of the LSDA 
and SIC. 
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II. LSDA SOLUTION 
A. Minimizing the energy functional 

We first solve the spin-polarized Kohn-Sham (KS) 
equations within the LSDA^ for an external potential 
set up by a Hydrogen ion and a sphere of positive charge 
centered on the ion with charge density uq and radius 
Rjeii- The system is overall charge neutral. 

We spherically symmetrize the spin densities, n'^(r), 
after each iteration of the self-consistency cycle before 
recalculating the KS potentials. The KS potential will 
therefore also be spherically symmetric and so the KS 
equations reduce to the radial Schrodinger equations. 
The self-consistent solutions^^ to these equations mini- 
mize the LSDA energy functional 

71,1,171, a 

1 f fn{r)n{r') ^ , , f , , , 

2 J J \r - r'\ ^ J 

Jr=o org 

where V'nZmli") ^'^s the self-consistent KS orbitals, vf{r) 
are the spherically symmetrized spin densities, n{r) ~ 
^^n°'(r) is the total electron density and Vext{r) is the 
external potential 

Vext[r) = / -dr 

The final two terms in the energy functional include 
the Coulomb repulsion between the ion and the positive 
background and the self-Coulomb repulsion of the posi- 
tive background. The electron density parameter, is 
defined by |7rr^no — 1 and N — |7ri?|g;;no is the number 
of electrons in the jellium sphere. 

The kinetic energy can be calculated either by numer- 
ically evaluating the radial derivatives of the orbitals, or 
alternatively by rearranging the KS equations and sub- 
stituting in for the Laplacian. In our calculations both 
methods are used and are found to give the same results. 

The energy functional is also used to calculate the en- 
ergy of the separate jellium sphere and Hydrogen atom. 
The immersion energy is then obtained by subtracting 
these quantities from the Hydrogen atom in jellium en- 
ergy. 

B. Jellium sphere sizes 

The number of electrons in the sphere is chosen by 
studying, within the LSDA, the immersion energy of 



a Hydrogen atom in a jellium sphere of fixed no as a 
function of the number of electrons in the sphere. Fig- 
ure [T] illustrates the dependence of the immersion en- 
ergy on the size of the jellium sphere for a background 
density of O.OOTa^'^ (r^ — 3.25). The most obvious fea- 
ture of the plot is the small scale bunching of immer- 
sion energies, with the immersion energy increasing or 
decreasing slightly as a particular angular momentum 
shell is filled. A larger scale feature is that the immer- 
sion energy oscillates around the immersion energy of 
Hydrogen in infinite jelliumjii^ These are Friedel oscil- 
lations and the wavelength predicted by the theory is 
AR/rg — n/{rskF) = 7r/(rs(37r^no)3 ) = 1.637 which is 
in good agreement with the wavelength as read off from 
Fig.m 

The amplitude of the Friedel oscillations in Fig. [T] be- 
comes smaller as the size of the jellium sphere is in- 
creased. This tells us that we can make the immersion 
energy arbitrarily close to the infinite jellium immersion 
energy by making the jellium sphere very large. However 
the damping of the oscillations is slow (unlike calcula- 
tions by Hintermann et al.^'^). Therefore because we are 
limited in the number of electrons we can include in the 
VQMC calculation we must exploit the cross-over points 
at which the immersion energy is as close as possible to 
the infinite jellium case. 

We observe that plots of immersion energy against 
Rjeii I , such as Fig. [H look the same for different choices 
of no. In particular the same sinusoidal oscillation is ob- 
served and the crossing points occur at the same values 
of Rjeii/rg — Ns, Therefore if for a given choice of no, 
one chooses a value of A^ for which the immersion energy 
is close to that of the infinite jellium, then one is assured 
that the immersion energy for all other values of no will 
also be close to the infinite jellium immersion energy. 

This is demonstrated by picking two values of A^ which, 
from Fig. [Tl have immersion energies close to the infinite 
jellium immersion energy. These values are 10 and 50, 
which are highlighted in Fig. [TJ In Fig. [5] the immersion 
energy is plotted for these two choices of A^ for a range 
of no values. The curve for the 10-electron jellium sphere 
is seen to give a reasonable approximation of the infinite 
jellium curve. However the 50-electron curve improves 
further by giving the minimum in the same place as for 
the infinite jellium. One sees that this curve is rigidly 
shifted above that of the infinite jellium, confirming that 
the error in the immersion energy due to the finite size 
is approximately independent of uq. 

Our calculations also show that the atom-induced den- 
sity of the Hydrogen atom in finite jellium tends smoothly 
towards that of a Hydrogen atom in infinite jellium as the 
size of the jellium sphere is increased. Figurc[2]shows our 
calculations of the atom-induced densities for a Hydro- 
gen atom in jellium spheres with a background density 
O.Ola^^. We consider jelfium spheres with 10, 50 and 338 
electrons and also plot the atom-induced density for a 
Hydrogen atom in infinite jellium. We see clearly that as 
we increase the number of electrons in the jellium sphere. 
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FIG. 1: The immersion energy versus jellium sphere radius for 
jellium spheres of background density 0.007o^^. The straight 
line is the value of the immersion energy for Hydrogen in infi- 
nite jellium of the same background density and the sinusoidal 
line is a guide to the eye. The highlighted points correspond 
to iV = 10 and N = 50. 
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FIG. 2: Plots of immersion energy versus background density 
for a Hydrogen atom immersed in jellium spheres of size 10 
(solid curve) and 50 (dashed curve) electrons. Also plotted is 
the immersion energy curve for a Hydrogen atom in infinite 
jellium (dot-dash curve). 



the atom-induced density approaches the limiting atom- 
induced density of a Hydrogen atom in infinite jellium. 



III. VQMC SOLUTION 
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FIG. 3: Plots of LSDA atom-induced densities for Hydrogen 
in finite jellium spheres of background density O.Olag"^, with 
(a) 10, (b) 50 and (c) 338 electrons (solid lines). Also plotted 
is the atom-induced density for a Hydrogen atom in infinite 
jellium at the same background density (dashed lines). The 
jellium sphere radius is shown by a vertical straight line. 



A. Methodology 



In VQMC we use importance sampling to replace the 
expectation value of the Hamiltonian with a sum over 
so-called local energies 



^'5.(R)i7«'T(R)dR 
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where M = j |^'T(R)P<iR, Ri are sets of particle coordi- 
r„(*) ^ . . . ^ r^'' } (referred to as configurations) 



nates, {rj*\ 

which sample the probability distribution l\l/r(R)P/A/' 
and where -ff\l/T(Ri)/^r(Ri) is referred to as the local 
energy, which is summed over a suitably large number 
of configurations, n. These configurations are generated 
using the Metropolis algorithmii^ 

By the variational principle we know that this expec- 
tation value is an upper bound on the exact ground-state 
energy. To get as close as possible to the ground-state, 
we minimize ai,e. ^^^^ the standard deviation of the local 
energy. 

The trial wavefunction, is written as a product of 
spin-up and spin-down Slater determinants, and , 
containing LSDA orbitals calculated for the same system. 
In addition the wavefunction contains a Jastrow factor— 
which describes the correlations between the electrons. 
For our Jastrow factor we use the same as that used 
by Sottile et alJ^ in their calculations of jellium spheres 
(they did not include an embedded atom), except for a 
'multipolar' term which we do not include in our calcu- 
lations. The trial wavefunction is 



^'T(ri, • • • ,rAr) = exp 




X exp 



T - 

\<i<j<N 



B. Fixing the jellium sphere size 

The number of electrons in the jellium sphere is chosen 
so that the immersion energy is as close to that of the 
atom in infinite jellium as possible. In addition, in order 
to have a good trial wavefunction we require a number 
of electrons such that it is possible to choose a configu- 
ration of the electrons in the open-shell so as to give a 
purely real wavefunction overall. The VQMC energy is 
also minimized by choosing configurations which are as 
non-magnetic as possible. With these considerations in 
mind, and also with a view to minimize computational 
effort, we choose a 10-electron jellium sphere for the cal- 
culations presented below. This choice leads to the fol- 
lowing VQMC electron configuration for the Hydrogen in 
jellium 

and the configuration for the 10-electron jellium sphere 
without the Hydrogen atom 

ls^\ls^^ ,2p^\2p^^ ,3d\3d^ 

where the d electrons are placed in orbitals with magnetic 
quantum number, m — 0. 

In the LSDA solution for the 10-electron jellium 
sphere, the two 3d electrons are instead filled as 3d^^ 
in order to minimize the total energy, in a manner con- 
sistent with Hund's rule. Note that the d sub-shell is not 
full and so the density contribution from these electrons 
is spherically symmetrized in LSDA before recalculating 
the KS potentials. 



IV. SIC SOLUTION 



1 + Cij {rj ) rij + dij r-ij ^ J _ 

xDT(ri, • • • ,rjv/2)i5'''(rjv/2+i, • • • ,rAr) 
where Cij is given by 

= cl + 4 arctan[(r2 - i?|,„)/2 A R^m] 

Parameters , bij , c^^ , and dij depend only on the 

relative spins of electrons i and j and depend only 
on the spin of electron i. The parameters aij are deter- 
mined by the electron-electron cusp condition. Note that 
the nuclear-electron cusp condition is automatically satis- 
fied on account of the LSDA orbitals in the determinant. 
The other 15 parameters in the Jastrow factor are varia- 
tional parameters which are varied in order to minimize 
CTj.e.. Correlated sampling'^'— i^l is used in the minimiza- 
tion procedure, which avoids the need to recalculate new 
configurations for each set of Jastrow parameters. 



In our SIC calculations, the SIC will not be applied to 
all electrons in the system. Instead we apply SIC only to 
the Is spin-up and spin-down electrons. This approach is 
taken because our Hydrogen atom in a jellium sphere is 
meant to be an approximation of Hydrogen in infinite jel- 
lium, and in the latter, only the Is^ and Is^ electrons are 
sufficiently localized to warrant the correction. Hence for 
our SIC solution, the Is electrons will obey KS equations 
with a Sl-corrected potential and all the other electrons 
obey the standard LSDA KS equations. 

For our 10-electron jellium sphere, the Is and 2s or- 
bitals of a given spin are calculated from different eigen- 
value equations (unlike in the LSDA, where both share 
the same potential) and so will not in general be orthog- 
onal. However, the KS formulation of DFT requires this 
orthogonality and so it must now be imposed, for which 
we use Gram- Schmidt orthogonalization^^ 

V'2r'M - AA- (^r2s{r) - / r,:{r')r2s{r')dv'r,s{r)^ 

where J^f^ ensures the correct normalization of the or- 
bitals. 
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The SIC energy functional is 



V- / ^-'■*''*(r)(--V2)^2T'^*''(r)dr+ 



n{r)n{r') 



drdr' 



E 



Zn^"!^dr^'-^ 
Jr=o r 5 

Note that the 2s contribution to the kinetic energy is 
most easily calculated by evaluating the radial derivatives 
of the orbitals. The alternative method of substituting 
in from the KS equations is also possible, but is more 
complicated in this case as the orthogonalized orbitals 
are now mixtures of the eigenstates of the KS equations. 
As with the LSDA results, both methods are found to 
give the same results. 



V. RESULTS 
A. Total energies and immersion energies 
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FIG. 4: Total energy of a Hydrogen atom immersed in a 10- 
electron jellium sphere for different background densities. 
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FIG. 5: Total energy of a lO-electron jellium sphere for dif- 
ferent background densities. 



We report VQMC solutions for a Hydrogen atom im- 
mersed in lO-electron jellium spheres of background den- 
sities O.OOIa^'^ through to O.OSa^'^ (r^ values of 6.2 and 2 
respectively). Total energies were calculated and are pre- 
sented, along with LSDA and SIC results, in Fig. 21 Total 
energies are also calculated using a close approximation 
to the Hartree-Fock method (HF). These points are cal- 
culated by evaluating the expectation value of the Hamil- 
tonian for a wavefunction consisting of just the Slater 
determinant part containing the LSDA orbitals. This is 
achieved by performing a VQMC run with the Jastrow 
factor set equal to one. 

Energies for jellium spheres but without the Hydrogen 
atom are presented in Fig. [5] Finally, the immersion 
energies are reported in Fig. [SI 

In calculating the immersion energy for HF, VQMC 
and SIC the exact value of the Hydrogen atom energy is 
used, namely — I3.606el^. This is appropriate as calcula- 
tions using HF and SIC both give this value and in the 
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FIG. 6: Immersion energies for a Hydrogen atom immersed in 
a 10-electron jellium sphere for different background densities. 
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case of VQMC, the calculation is correct to a very high 
level of precision. For the LSDA calculation, the atom 
energy is found to be — 13.030eV. 

The total energy curves calculated using the different 
methods are all of the same shape and feature minima at 
roughly the same locations. The LSDA, SIC and VQMC 
curves lie within leV of one another in the atom in jel- 
lium total energy graph and similarly the LSDA and 
VQMC curves are within leV of one another in the jel- 
lium sphere total energy graph. HF gives total energy 
curves for the atom in jellium and for the jellium sphere 
which are rigidly shifted by about 5cV above the VQMC 
curves. 

We note that our VQMC total energy curve for the 
jellium sphere is slightly different to the type obtained 
by Sottile et alM Their VQMC results were found to be 
slightly lower than the LSDA results above background 
densities of about O.Ola^'^. This is because we neglected 
to include the 'multipolar' term in the trial wavefunction 
which Sottile et alr^ used in their calculations. In their 
work, this term was found to improve the VQMC solution 
for the larger background densities, but was relatively 
unimportant at densities near the minimum. 

When the total energy curves are subtracted from one 
another to give the immersion energy curves, the large 
error in the HF calculations cancel out. We find all four 
curves lie within leV of one another and all curves feature 
a minimum at approximately 0.004a^^. Notice that the 
LSDA results for the immersion energy curve are lower 
than the VQMC results, which means that the LSDA is 
overbinding the atom to the jellium relative to VQMC. 

Calculations by Sugiyama et alJ^ of the immersion en- 
ergy has a discrepancy of 5eV between Monte Carlo and 
LSDA, which is very much larger than we find here. Pos- 
sibly this difference was because of a different choice of 
trial wavefunction. 

We see that the VQMC immersion energy curve below 
0.005a^^ is more closely followed by the SIC curve than 
by the LSDA curve. If we regard our VQMC results as a 
benchmark (see Section I VII for a discussion of this) then 
this indicates that the SIC is outperforming the LSDA 
for these background densities. This is expected as the 
Is bound states are becoming more localized at these low 
background densities and therefore the self-interaction of 
these bound states is increasing. 

Note that for an atom in infinite jellium the SIC and 
LSDA curves would begin to coincide at large background 
densities due to the Is electrons becoming increasingly 
delocalized. However for our system of an atom in finite 
jellium this is not the case. As the background density is 
increased, the sphere size becomes smaller and so effects 
due to the finite size of our system begin to impinge on 
the results. 
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FIG. 7: Atom-induced density as calculated using the differ- 
ent methods, for no — 0.004ag"^. The solid line is the LSDA 
density, the dashed line is the SIC density and the noisy line 
is the VQMC density. 



B. Atom-induced densities 

Atom-induced densities (multiplied by r^) are plotted 
in Fig. [7| for a background density of 0.004a^^. The 
Hydrogen-like maximum in the atom-induced density is 
very similar between VQMC and LSDA for r less than 
AaB- A discrepancy of at most 5% is identified in this re- 
gion. On the other hand, as highlighted by Hoffman and 
Pratt the work of Sugiyama et alJ^ shows a difference 
of 30% between LDA and VQMC atom-induced densities 
near the proton. The much smaller discrepancy in our 
results is probably due to our use of a trial wavefunction 
which performs better than that of Sugiyama et al. close 
to the proton. 

However, our atom-induced density plots show the 
same discrepancy in the Friedel oscillation between LSDA 
and VQMC as was seen by Sugiyama et al ^^^'^^ In par- 
ticular the minimum in the atom-induced density for the 
LSDA is much deeper than that for the VQMC. Also the 
wavelength of the oscillation is larger for VQMC and so 
the Friedel oscillation maximum occurs at a larger radius. 

The SIC atom-induced density is also included in 
Fig. [71 The SIC and LSDA densities are essentially iden- 
tical, which seems unusual given the difference in the 
immersion energies for these methods. 



VI. CONCLUSIONS 

We have presented what we believe to be the first re- 
liable VQMC results for the atom in jellium model sys- 
tem. The immersion energy versus background density 
curves for VQMC, LSDA and SIC lie within leV of each 
other with minimas in approximately the same positions. 
The LSDA results show an overbinding of the atom to 
the jellium relative to the VQMC result. Viewing the 
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VQMC as a benchmark, this is consistent with the gen- 
eral overbinding seen in LSDA. In addition, for low back- 
ground densities, the immersion energy curve of VQMC 
is more closely followed by the the SIC immersion energy 
curve than by the LSDA immersion energy curve. Again, 
viewing VQMC as a benchmark, this is consistent with 
the fact that the SIC is expected to be more accurate 
than the LSDA for systems with more strongly localized 
electrons (as is the case for the low background densities). 

The status of the VQMC results as a benchmark is 
uncertain however. More accurate DQMC calculations 
would yield lower total energies of both the atom in jel- 
lium and the jellium relative to the VQMC results. How- 
ever, because the immersion energy is calculated as the 
difference between these energies, part of the change in 



going from VQMC to DQMC will cancel out when we 
calculate the immersion energy. Furthermore, Sottile et 
alM found a rigid shift of only -0.2eV in the DQMC 
energies of an 8-electron jellium sphere relative to the 
VQMC energies. This suggests that our results for the 
immersion energy would not be strongly modified by us- 
ing DQMC. 
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